/* Copyright 2007 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 * Marcus Cuda
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;
using dnAnalytics.Properties;

namespace dnAnalytics.Math
{
    /// <summary>
    /// This class implements a collection of special function evaluations. This class has a static constructor which 
    /// will precompute a small number of values for faster runtime computations.
    /// </summary>
    public static partial class SpecialFunctions
    {
        // Constants used in the computation of the gamma function.
        private const int Gamma_n = 10;
        private const double Gamma_r = 10.900511;
        private static readonly double Gamma_cst = 2.0 * System.Math.Sqrt(System.Math.E / System.Math.PI);

        private static readonly double[] Gamma_dk = new double[]
                                                    {
                                                        2.48574089138753565546e-5,
                                                        1.05142378581721974210,
                                                        -3.45687097222016235469,
                                                        4.51227709466894823700,
                                                        -2.98285225323576655721,
                                                        1.05639711577126713077,
                                                        -1.95428773191645869583e-1,
                                                        1.70970543404441224307e-2,
                                                        -5.71926117404305781283e-4,
                                                        4.63399473359905636708e-6,
                                                        -2.71994908488607703910e-9
                                                    };

        private static readonly double GammaLn_cst = System.Math.Log(2.0) + 0.5 * System.Math.Log(System.Math.E / System.Math.PI);
        private static readonly double GammaLn_logpi = System.Math.Log(System.Math.PI);

        
        /// <summary>
        /// Static constructor for the special functions class. This method precomputes a small number of values for faster runtime computation.
        /// </summary>
        static SpecialFunctions()
        {
            // Precompute factorial values that are representable in double format.
            Factorial_cache = new double[Factorial_max + 1];
            Factorial_cache[0] = 1.0;
            for (int i = 1; i <= Factorial_max; i++)
            {
                Factorial_cache[i] = Factorial_cache[i - 1] * i;
            }
        }

        /// <summary>
        /// Computes the logarithm of the Gamma function. 
        /// </summary>
        /// <param name="z">The argument of the gamma function.</param>
        /// <returns>The logarithm of the gamma function.</returns>
        /// <remarks>
        /// This implementation of the computation of the gamma and logarithm of the gamma function follows the derivation in
        ///     "An Analysis Of The Lanczos Gamma Approximation", Glendon Ralph Pugh, 2004.
        /// We use the implementation listed on p. 116 which achieves an accuracy of 16 floating point digits. Although 16 digit accuracy
        /// should be sufficient for double values, improving accuracy is possible (see p. 126 in Pugh).
        /// </remarks>
        public static double GammaLn(double z)
        {
            if (z < 0.5)
            {
                double s = Gamma_dk[0];
                for (int i = 1; i <= Gamma_n; i++)
                {
                    s += Gamma_dk[i] / (i - z);
                }
                return GammaLn_logpi - System.Math.Log(System.Math.Sin(System.Math.PI * z)) - System.Math.Log(s) - GammaLn_cst - (0.5 - z) * System.Math.Log((0.5 - z + Gamma_r) / System.Math.E);
            }
            else
            {
                double s = Gamma_dk[0];
                for (int i = 1; i <= Gamma_n; i++)
                {
                    s += Gamma_dk[i] / (z + i - 1.0);
                }
                return System.Math.Log(s) + GammaLn_cst + (z - 0.5) * System.Math.Log((z - 0.5 + Gamma_r) / System.Math.E);
            }
        }

        /// <summary>
        /// Computes the Gamma function. 
        /// </summary>
        /// <param name="z">The argument of the gamma function.</param>
        /// <returns>The logarithm of the gamma function.</returns>
        /// <remarks>
        /// This implementation of the computation of the gamma and logarithm of the gamma function follows the derivation in
        ///     "An Analysis Of The Lanczos Gamma Approximation", Glendon Ralph Pugh, 2004.
        /// We use the implementation listed on p. 116 which achieves an accuracy of 16 floating point digits. Although 16 digit accuracy
        /// should be sufficient for double values, improving accuracy is possible (see p. 126 in Pugh).
        /// </remarks>
        public static double Gamma(double z)
        {
            if (z < 0.5)
            {
                double s = Gamma_dk[0];
                for (int i = 1; i <= Gamma_n; i++)
                {
                    s += Gamma_dk[i] / (i - z);
                }
                return System.Math.PI / (System.Math.Sin(System.Math.PI * z) * s * Gamma_cst * System.Math.Pow((0.5 - z + Gamma_r) / System.Math.E, 0.5 - z));
            }
            else
            {
                double s = Gamma_dk[0];
                for (int i = 1; i <= Gamma_n; i++)
                {
                    s += Gamma_dk[i] / (z + i - 1.0);
                }
                return s * Gamma_cst * System.Math.Pow((z - 0.5 + Gamma_r) / System.Math.E, z - 0.5);
            }
        }

        /// <summary>
        /// Computes the logarithm of the Euler Beta function.
        /// </summary>
        /// <param name="z">A positive real number.</param>
        /// <param name="w">A positive real number.</param>
        /// <returns>The logarithm of the Euler Beta function evaluated at z,w.</returns>
        /// <exception cref="ArgumentException">If <paramref name="z"/> or <paramref name="w"/> are not positive.</exception>
        public static double BetaLn(double z, double w)
        {
            if (z <= 0.0)
            {
                throw new ArgumentException(Resources.NotPositive, "z");
            }
            if (w <= 0.0)
            {
                throw new ArgumentException(Resources.NotPositive, "w");
            }
            
            return GammaLn(z) + GammaLn(w) - GammaLn(z + w);
        }

        /// <summary>
        /// Computes the Euler Beta function.
        /// </summary>
        /// <param name="z">A positive real number.</param>
        /// <param name="w">A positive real number.</param>
        /// <returns>The Euler Beta function evaluated at z,w.</returns>
        /// <exception cref="ArgumentException">If <paramref name="z"/> or <paramref name="w"/> are not positive.</exception>
        public static double Beta(double z, double w)
        {
            return System.Math.Exp(BetaLn(z, w));
        }

        // This is the largest integer for which the factorial function in double format doesn't overflow.
        private const int Factorial_max = 170;
        static private readonly double[] Factorial_cache;

        /// <summary>
        /// Computes the factorial function: x!. The function can represent all number up to 22! exactly, all numbers
        /// up to 170! using a double representation. All larger values will overflow.
        /// </summary>
        /// <param name="x">A nonnegative integer for which to compute the factorial function. </param>
        /// <returns>The factorial of x: x!</returns>
        /// <exception cref="ArgumentException">If <paramref name="x"/> is negative.</exception>
        public static double Factorial(int x)
        {
            if (x < 0)
            {
                throw new ArgumentException(Resources.NotPositive, "x");
            }
            if (x <= Factorial_max)
            {
                return Factorial_cache[x];
            }
            return Double.NaN;
        }

        /// <summary>
        /// Returns the log of the factorial function.
        /// </summary>
        /// <param name="x">A nonnegative integer for which to compute the logarithm of the factorial function.</param>
        /// <returns>The logarithm of the factorial function.</returns>
        /// <exception cref="ArgumentException">If <paramref name="x"/> is negative.</exception>
        static public double FactorialLn(int x)
        {
            if (x < 0)
            {
                throw new ArgumentException(Resources.NotPositive, "x");
            }
            if (x <= 1)
            {
                return 0.0;
            }
            return GammaLn(x + 1.0);

        }

        /// <summary>
        /// Computes the digamma function which is mathematically defined as the derivative of the logarithm of the gamma function.
        /// This implementation is based on
        ///     Jose Bernardo
        ///     Algorithm AS 103:
        ///     Psi ( Digamma ) Function,
        ///     Applied Statistics,
        ///     Volume 25, Number 3, 1976, pages 315-317.
        /// Using the modifications as in Tom Minka's lightspeed toolbox.
        /// </summary>
        /// <param name="x">The argument of the digamma function.</param>
        /// <returns></returns>
        static public double DiGamma(double x)
        {
            const double c = 12.0,
                d1 = -0.57721566490153286,
                d2 = 1.6449340668482264365,
                s = 1e-6,
                s3 = 1.0/12.0,
                s4 = 1.0/120.0,
                s5 = 1.0/252.0,
                s6 = 1.0/240.0,
                s7 = 1.0/132.0;

            if ( System.Double.IsNegativeInfinity(x) || System.Double.IsNaN(x) )
            {
                return System.Double.NaN;
            }

            // Handle special cases.
            if ( x <= 0 && System.Math.Floor(x) == x )
            {
                return System.Double.NegativeInfinity;
            }

            // Use inversion formula for negative numbers.
            if (x < 0)
            {
                return DiGamma(1.0 - x) + System.Math.PI / System.Math.Tan(-System.Math.PI * x);
            }

            if(x <= s)
            {
                return d1 - 1/x + d2*x;
            }

            double result = 0;
            while(x < c) {
                result -= 1/x;
                x++;
            }

            if(x >= c)
            {
                double r = 1/x;
                result += System.Math.Log(x) - 0.5*r;
                r *= r;

                result -= r * (s3 - r * (s4 - r * (s5 - r * (s6 - r * s7))));
            }

            return result;
        }

        /// <summary>
        /// Computes the binomial coefficient: n choose h.
        /// </summary>
        /// <param name="n">A nonnegative value n.</param>
        /// <param name="h">A nonnegative value h.</param>
        /// <returns>The binomial coefficient: n choose h.</returns>
        /// <exception cref="ArgumentException">If <paramref name="n"/> or <paramref name="h"/> are negative.</exception>
        public static double Binomial(int n, int h)
        {
            if (n < 0)
            {
                throw new ArgumentException(Resources.NotPositive, "n");
            }
            if (h < 0)
            {
                throw new ArgumentException(Resources.NotPositive, "h");
            }
            return System.Math.Floor(0.5 + System.Math.Exp(FactorialLn(n) - FactorialLn(h) - FactorialLn(n - h)));
        }

        /// <summary>
        /// Computes the logarithm of the binomial coefficient: n choose h.
        /// </summary>
        /// <param name="n">A nonnegative value n.</param>
        /// <param name="h">A nonnegative value h.</param>
        /// <returns>The log of the binomial coefficient: n choose h.</returns>
        /// <exception cref="ArgumentException">If <paramref name="n"/> or <paramref name="h"/> are negative.</exception>
        public static double BinomialLn(int n, int h)
        {
            if (n < 0)
            {
                throw new ArgumentException(Resources.NotPositive, "n");
            }
            if (h < 0)
            {
                throw new ArgumentException(Resources.NotPositive, "h");
            }
            return FactorialLn(n) - FactorialLn(h) - FactorialLn(n - h);
        }

        /// <summary>
        /// Computes the multinomial coefficient: n choose n1, n2, n3, ...
        /// </summary>
        /// <param name="n">A nonnegative value n.</param>
        /// <param name="ni">An array of nonnegative values that sum to <paramref name="n"/>.</param>
        /// <returns>The multinomial coefficient.</returns>
        /// <exception cref="ArgumentNullException">if <paramref name="ni"/> is <see langword="null" />.</exception>   
        /// <exception cref="ArgumentException">If <paramref name="n"/> or any of the <paramref name="ni"/> are negative.</exception>
        /// <exception cref="ArgumentException">If the sum of all <paramref name="ni"/> is not equal to <paramref name="n"/>.</exception>
        public static double Multinomial(int n, int[] ni)
        {
            if (n < 0)
            {
                throw new ArgumentException(Resources.NotPositive, "n");
            }
            if (ni == null)
            {
                throw new ArgumentNullException("ni");
            }
 
            int sum = 0;
            double ret = FactorialLn(n);
            for (int i = 0; i < ni.Length; i++)
            {
                if (ni[i] < 0)
                {
                    throw new ArgumentException(Resources.NotPositive, "ni["+i+"]");
                }

                ret -= FactorialLn(ni[i]);
                sum += ni[i];
            }
            
            // Before returning, check that the sum of all elements was equal to n.
            if (sum != n)
            {
                throw new ArgumentException(Resources.ParameterNotConformable, "ni");
            }

            return System.Math.Floor(0.5 + System.Math.Exp(ret));
        }

        /// <summary>
        /// Computes the <paramref name="t"/>'th Harmonic number.
        /// </summary>
        public static double Harmonic(int t)
        {
            return Constants.EulerMascheroni + DiGamma(t + 1.0);
        }
    }
}
